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We present an algorithm to simulate the many-body depletion interaction between anisotropic colloids in an 
implicit way, integrating out the degrees of freedom of the depletants, which we treat as an ideal gas. Because 
the depletant particles are statistically independent and the depletion interaction is short-ranged, depletants 
are randomly inserted in parallel into the excluded volume surrounding a single translated and/or rotated 
colloid. A configurational bias scheme is used to enhance the acceptance rate. The method is validated and 
benchmarked both on multi-core CPUs and graphics processing units (GPUs) for the case of hard spheres, 
hemispheres and discoids. With depletants, we report novel cluster phases, in which hemispheres first assemble 
into spheres, which then form ordered hcp/fcc lattices. The method is significantly faster than any method 
without cluster moves and that tracks depletants explicitly, for systems of colloid packing fraction < 0.50, 
and additionally enables simulation of the fluid-solid transition. 


I. INTRODUCTION 

The self-assembly of anisotropic particles into complex 
structures has emerged as a promising strategy towards 
the fabrication of materials with novel propertied Meth¬ 
ods for the synthesis of anisotropic nano- and colloidal 
particled^ are becoming available, and enable experi¬ 
ments that study their phase behavioi®^. Anisotropic 
particles, such as proteins, are also emerging build¬ 
ing blocks for biomateriald^. Simulations predict a 
wealth of different crystal structures t hat hard shapes 
form through maximization of entropy.^^E^ In addition 
to particle shape, attractive interactions between patchy 
particles can be important in achieving desired target 
structured^EIl, Towards that end, the main routes that 
are actively being explored include surface funct ional- 
ization of nanoparticles using short DNA moleculeJ^^^i^, 
and exploiting the depletion interaction betwee n colloids 
in the presence of small polymer chaindE^E^. Here we 
focus on the depletion interaction, since it is of entropic 
origin and arises without the need for engineering par¬ 
ticle surface chemistry, emerging in mixtures of colloids 
with non-adsorbing polymer. 

DepletiojJi^ describes the emergent attraction between 
colloids in solution that maximize the free volume avail¬ 
able to a small-particle cosolute via overlap of their ex¬ 
cluded volume shells. It has been demonstrate d th at de¬ 
pletion enhances the directional entropic force^^lMMl j.g_ 
suiting from anisotropic particle shape, and that it pro¬ 
motes the contact between large facets. The depletion 
interac tion can promote binding between lock and key 
colloid^^l2U and lead to the formation of porous phase^^. 
Because depletion mediates an additional attraction of 
entropic origin, this interaction can be thought of as 
competing with contact (excluded volume) interactions 
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FIG. 1. Explicit {left) vs. implicit (right) treatment of de¬ 
pletion interactions. Hard tetrahedra in solution with small, 
penetrable hard spheres aggregate face to face, to maximize 
the free volume available to the depletants. 


resulting from particle shape. Depletion thus enables 
novel phase behavior through the additional parameters 
of depletant shape and densitj®^. Therefore, it is desir¬ 
able to have a method to investigate the self-assembly of 
anisotropic shapes in the presence of depletants. Results 
for the phase behavior of binary hard sphere mixtures 
have been reportecPSl using thermodynamic integration. 
In general, however, such results are challenging to ob¬ 
tain because of the size disparity between the colloid and 
the depletant. If one is interested in the phase behavior 
of the colloids, a customary ap proxi mation treats the de¬ 
pletant particles as an ideal ga d^^ l ^^ l This approximation 
would, in principle, allow integrating out the depletant to 
arrive at an effective colloid-colloid interaction; however, 
the resulting interaction is a many-body interaction and 
we are not aware of any prior implementation that treats 
many-body effects exactly. Here, we propose a novel, 
parallel Monte Carlo algorithm to simulate the depletion 
interaction between arbitrarily shaped colloids in an efh- 
cient manner that includes many-body effects. 

Figure shows the effect of depletion interactions be¬ 
tween two hard tetrahedra in solution with small pen¬ 
etrable hard spheres. The small spheres mediate an at- 
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tractive interaction between the colloids that drives them 
to aggregate face to face. For two particles only, the de¬ 
pletion interaction can be easily simulated explicitly (left 
panel) or implicitly (right panel). However, implicit sim¬ 
ulation of depletion interactions allow for a tremendous 
performance benefit, particularly for dilute systems of 
colloids and high densities of depletants, as we demon¬ 
strate below. 

This paper is organized as follows. In section |TTj we 
discuss previous numerical methods for the simulation 
of depletion interactions. We describe our algorithm in 
section [IlIl and validate it against published data for hard 
spheres in the following section HYl Section [V] contains 
new results for hemispheres and discoid^^in the presence 
of depletants, obtained with the new algorithm. Finally, 
in Sec. |VI| we summarize and give an outlook on future 
applications of the method. 


II. BACKGROUND 

Previous numerical treatments of depletion interac¬ 
tions employ cluster moves. Biben, Bol huis a nd Frenkel 
proposed a configurational bias approac H^^F8 ( where de¬ 
pletants overlapping with a moved colloid are reinserted 
to enhance the acceptance probability of colloid moves. 
A geometric cluster algorithm has also been proposed 
by Dress and KrauttP^, which is rejection-free and can 
therefore greatly enhance the equilibration of dilute sys¬ 
tems of colloids. However, when the system is dense in 
colloids, clusters can ^an the system and the algorithm 
ceases to be efficientPSl. To explore the phase behavior 
of a system of hard spheres in penetrable hard-sphere 
depletants, Vink and Horbach proposed grand-canonical 
simulation of both the colloids and the depletants, and 
they could efficiently sample the gas-liquid coexistence 
curv^^. However, their scheme does not generalize well 
beyond to the fluid-solid transition, because it is based 
on particle insertion. 

All these methods have in common that they track the 
small depletant particles explicitly, which are stored in 
memory. An interesting alternative was proposed by Di- 
jkstra et al.^^, who proposed a Monte Carlo integration 
of the free volume around every single moved colloid. 
However, their scheme does not obey detailed balance, 
and achieving sufficient accuracy comes at the expense 
of computation time, as we discuss in more detail be¬ 
low. Another implicit implementation of the depletion 
interaction between octahedra was proposed by Henzie 
et al.l^, where the generally anisotropic many-body in¬ 
teraction is reduced to an isotropic pair potential. We 
note that such a drastic simplification, while rendering 
the problem computationally tractable, is insufficient to 
allow the study of arbitrary shapes. 

The scheme we describe in the following section is a 
completely general treatment of depletion interactions 
between anisotropic particles due to an ideal gas of deple¬ 
tants, and works well both for dilute and dense systems. 


In the ideal gas treatment, depletants interact with col¬ 
loids but not with each other. The algorithm is rigorous, 

i.e. it obeys detailed balance, and it can be efficiently 
implemented on multi-core processors and graphics pro¬ 
cessing units (GPUs). 


III. DESCRIPTION OF THE ALGORITHM 
A. Semigrand N/XpVT ensemble 

We simulate a semigrand ensemble of N colloids in a 
grand-canonical bath of penetrable depletants of chem¬ 
ical potential /ip. The partition sum for the depletants 
is 
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where Vf = Vf[fc,i] is the free volume available to de¬ 
pletants and Ap the thermal de Broglie wavelength as¬ 
sociated with the depletants. We denote the colloid- 
colloid contribution to the Hamiltonian as Hcc = 
Z)i,jGcolloids > where C/y = oo for two colloids that 
overlap, and Uij = 0 otherwise. The colloid-polymer 
contribution to the Hamiltonian Hep is defined analo¬ 
gously. Summation over the number Np of depletants in 
the system results in 

g-/3H{ry,,} = g^pV-/3ffcc^ (3) 


where Zp = 


is the depletant fugacity. 


B. Basic idea 

Our central algorithmic result is the following Monte 
Carlo scheme to integrate the colloids under the action 
of the effective potential Hes = —P~^ZpVf[fc,i] occurring 
in Eq. ([^. The basic idea of the algorithm, which we 
present here, is very simple, and we describe optimized 
versions of it in ensuing sections. 

1. Propose a trial move for the colloids M —>■ M'. 

2. Generate Np random depletant positions uni¬ 
formly in the free volume of the old configuration 
M, where Np is chosen according to PzpVfi^p) ~ 
Poisson(V/Zp), where Poisson(A) is the Poisson dis¬ 
tribution of mean and variance A. One possibil¬ 
ity is to use rejection sampling in a larger volume 
Uo D Vf. 

3. Reject the trial move if any depletant overlaps with 
the new colloid configuration M', otherwise accept. 
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FIG. 2. Depletant positions (disks) considered for rejection 
of a colloid move (shaded squares). The difference between 
configurations M and M' is the position of the dark shaded 
colloid. When moving the colloid to the new position, de- 
pletants are randomly inserted into the circumsphere of the 
excluded volume, and depletants that only overlap with the 
shape in the new configuration M' lead to rejection. Deple¬ 
tants that overlap with the colloid in the old position or with 
surrounding colloids (light shaded square) are not considered. 


In other words, we have an a priori move generation 
probability 

(M ^ M') = (M ^ {Np) (4) 

= P-UM ^ 

where P““(M —>■ M') is symmetric in Afyy O 
In Eq. Q, we have used the definition of the Poisson 
distribution Pzj,Vf{P^p) with average ZpVf, the number of 
depletants in the free volume. We impose the following 
acceptance probability 

M') = min(l, . (5) 


Figure contains a graphical summary of the algo¬ 
rithm. Here, a square colloid is moved from configuration 
M to configuration M', by some translation and/or rota¬ 
tion, and depletants are placed in the free volume. As we 
detail below in Sec. Ill C| the sampling can be restricted 
to the circle (or sphere, in three dimensions) containing 
the colloid in the new colloid position. By using rejec¬ 
tion sampling, any depletants falling into the excluded 
volume at the old position are ignored. Depletants that 
overlap only in the new configuration lead to a rejection 
of the colloid move. 

Next, we show that the above scheme obeys detailed 
balance, which is required for correctly sampling the en¬ 
semble defined by Eq. (|^ in the statistical sense. The 
transition probability tt from the old configuration M to 
the new configuration M' obeys 


^ M')Pi^P{M ^ M') 

lyp. 

min(I, ( 5 ) 


We require for detailed balance that ttm-j-m' = 
TTM'-i-M, and average over all realizations (^Np, 
of depletants, in the free volume V/, 


r dr^ 
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Note that in order to obtain Eq. Q, we observe that 
the Poisson distribution is normalized in such a way so as 
to cancel out the depletant contribution, to the en¬ 

semble weight. The integrand in the last line of Eq. Q 
is non-zero exactly for rp^. S P/; hence, after perform¬ 
ing the summation over Np, the transition probability 
becomes 


= e ^^“P///j(M ^ M')min(l,e 

^zpP(Vfnv})^ (g) 

where the volume /i(P/nP/) is the intersection of the free 
volume Vf in the old configuration and the free volume P/ 
in the new configuration. This term arises because of the 
integration domain in Eq. Q. Because of the symmetry 
of the Metropolis criterion, 


e ^^“‘=min(I,e = e ^^“'=min(e ^'^'^==,1) (9) 


and the symmetry property of the set intersection, the 
product in Eq. ([^ is symmetric under the exchange M 0 
M'. Consequently, our integration scheme obeys detailed 
balance. 


C. Improved formulation 


The above integration scheme conveys the general idea 
of the algorithm. However, this algorithm is impractical 
to implement as is in an actual program, because it would 
require computation of the free volume P/ in the entire 
simulation box for every single colloid move. Without 
loss of generality, we can restrict the sampling volume P/ 
for depletants to a smaller volume Vq D V/^ci \14xci, i-e. 
containing the excluded volume of the colloids in the 
system in the new configuration minus the excluded vol¬ 
ume Pexci in the old configuration. The improved scheme 
is the same as the old scheme (Sec. HIB), as are the 
move generation and acceptance probabilities, with the 
exception that Vf is replaced by P/ n Pq. The proof of 
detailed balance is only slightly more complicated for this 
algorithm. 

We rewrite the ensemble weight 


_ g-/9ifcc-l-ZpV/ 

^ ^-l3H^a+zp[p(VfnVo)+p{Vfn%)] 


( 10 ) 
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where Vq denotes the complement 14\Vo with respect 
to the simulation volume V. Using Eq. (10), integrat¬ 


ing over Vq n Vf and using transformations analogous to 
Eqs. ([^-([^, the transition probability M —>■ M' aver¬ 
aged over the number of test depletants and their posi¬ 
tions becomes 


= e-^^-min ^ M') 


^z^[tj.{VfnVonv})+tJ.iVfnVo)] _ 


( 11 ) 


It is straightforward to show that this transition prob¬ 
ability is symmetric for forward and reverse moves. Since 
Vq 3 I4x(.i\I4xci, it follows that 


K) C C^AKxcl C U Uexcl = V} U Uexcl (12) 

and therefore Vq = Vq n (Uj U I4xci)- Hence, applying the 
distributive law, 

VfnVo = VfnVon{v}u Kxci) = U/ni^nu/, (i3) 


because Vf n t4xci = 0- 

Using Eq. (131 we rewrite the transition probability 
Eq. (ITT) as 


= e 


min (1, M') 

gZp[fi{VfnvjnVo)+ti{VfnVfnVo)] 



overlapping at new position cannot be inserted 

(^ reinserted at old position ignored 

FIG. 3. Computation of the configurational bias weight for 
the forward move. When a single moved colloid overlaps with 
a randomly inserted depletant in the new configuration M', 
we attempt to reinsert it ritriai times such that it overlaps 
with the shape in the old conhguration M. Valid insertion 
attempts are those where the depletant neither overlaps with 
a surrounding colloid nor with the colloid in the new position. 
The configurational bias weight is computed from the number 
of successful reinsertions, cf. Eq. (|16[). 


excluded volume shell around the colloid particle on av¬ 
erage, moves will be rejected most of the time. Equi¬ 
libration of colloids in very dense depletant systems is 
therefore difficult. 


and because the measures in the exponent are taken from 
disjoint sets we can simplify this equation as 

= e-^^-min ^ M') 

e^pt^iVfnv}) ^2.5) 

This is the same transition rate as Eq. ([^, consequently 
our restricted sampling algorithm obeys detailed balance. 

We may choose Vq as the smallest region with Vq T 
VAci\Kixci that is convenient to sample from. E.g., we 
can sample in the excluded volume lAxd i of the sin¬ 
gle moved colloid i at the position of the new config¬ 
uration M' only, ignoring depletants that overlap with 
the colloid in the old configuration M. For anisotropic 
colloids, we will choose the circumsphere of diameter 
dcoiioid + ddepietant around the colloid in the new con¬ 
figuration, as done in Fig. 

We remark that a further possible optimization con¬ 
sists in restricting the sampling to the excluded volume 
shell of the moved colloid V^xci.AKore ^^.d it can be 
shown, using steps analogous to above, that such a choice 
also fulfills detailed balance. 


D. Configurational bias moves 

The algorithm described above gives finite acceptance 
rates for translation step sizes <5 A where R is 

the size of the colloid, which is in general anisotropic. 
However, when there is more than one depletant in the 


To ameliorate this situation, we apply the co nfigu ra- 
tional bias move of Biben, Bolhuis and Frenkep32Sl 
implicit depletants, the idea of which we briefly summa¬ 
rize. Figure [^depicts the basic idea. For every depletant 
overlapping in the new configuration M', we attempt to 
reinsert it ntriai times such that it overlaps with the shape 
in the old configuration M, but does not overlap with any 
other colloid. Such a cluster move obeys detailed balance 
because when performing the reverse move from M' to 
M, the reinserted colloid will overlap in the old configu¬ 
ration. To correct for the configurational bias generated 
in this wajEH, we modify the acceptance probability 


Hacc = min 


AT. 


overlap 


n 


(A^i„.e.t.i + l)iV, j ’ 


(16) 


in which insert,i and iVAggrt i number of times 

the overlapping depletant i could be reinserted without 
overlap into the old and new configuration, respectively. 
The numbers Ni,N[ < ntriai, count the valid insertion 
attempts in which the depletant overlaps with the moved 
shape in the old (new) configuration, without overlapping 
in the other. All other insertion attempts are ignored. 
The increment of one (iVinsert + 1) is necessary because 
the depletant the colloid was overlapping with originally 
counts as a successful reinsertion attempt for the reverse 
move. 
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E. Parallel implementation 

An important feature of our algorithm is that the de- 
pletant insertions are independent and can be performed 
in parallel. We exploit this feature to implement the 
algorithm on the GPU. Some details of the GPU imple¬ 
mentation are described in App. 

In addition, depletants are inserted only in a local 
neighborhood of the particle, reflecting the short-ranged 
nature of the depletion interaction. This means the 
parallelization scheme for particle based Monte Carlo 
that has recently been introduced wit hin th e Hard Par¬ 
ticle Monte Carlo (HPMC) frameworlPMl in HOOMD- 
biu^HHm can be generalized to our implicit depletion al¬ 
gorithm. HPMC uses a checkerboard decomposition to 
allow parallelization of the MC simulation on a graphics 
processor (GPU). The checkerboard is colored in such a 
way that simultaneously active cells are separated by a 
layer of inactive cells of width dcoiioid + ddepietantj which 
allows the active cells to be updated independently. Par¬ 
ticles are not allowed to move outside their cells. The 
checkerboard coloring is permuted randomly. In order to 
maintain ergodicity, the grid lines are randomly shifted. 
HPMC also runs on the CPU, using an efficient tree- 
based particle data storage for overlap checks in combi¬ 
nation with a sequential algorithm. Both the CPU and 
the GPU code path can be combined with spatial do¬ 
main decompositioiP^, using the same same concept of 
an inactive layer for parallel execution. A reference im¬ 
plementation of the algorithm described in this paper will 
be released open-source as part of HOOMD-blue^. 

IV. VALIDATION 

A. Equation of state of the penetrable hard sphere model 

To validate our method, we compare results for hard 
spheres with the previously obtained results by Dijkstra 
et al.^. We note that even though their implicit algo¬ 
rithm for depletion does not obey detailed balance, it 
relies on minimizing errors from the violation of detailed 
balance through increasing the discretization of the MC 
integration step, which is a trade-off between accuracy 
and performance. In order to obtain an accurate equa¬ 
tion of state, Dijkstra et. al had to restrict themselves to 
fairly small systems of N = 128 spheres. Fig. [^compares 
results obtained with our algorithm (filled symbols) to 
those from Fig. 2 of Ref. (stars). We show the mea¬ 
sured free volume fraction (fip available to the penetrable 
hard spheres of same size, as a function of the reservoir 
volume fraction (j)p for different colloid volume fractions 
(f>c at constant simulation volume. For a system size of 
N = 128 colloids, our and Dijkstra’s results are in es¬ 
sentially perfect agreement, mutually validating both al¬ 
gorithms (top panel). However with our new algorithm 
we can easily perform simulations for a larger system of 
N = 1000 spheres. We do see slight deviations from the 
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FIG. 4. Equation of state of spheres in penetrable hard sphere 
depletants. Plotted is the measured free volume 4>p available 
to the penetrable hard spheres of size ratio q = ddep/dcoiioid = 
1 vs. the reservoir volume fraction (pp of the depletants, for 
different hard sphere volume fractions = 0.01. ..0.3 (filled 
symbols). Data by Dijkstra et al.^^for N = 128 is shown as 
asterisks. Upper panel: equation of state for N — 128 colloids, 
lower panel: N = 1000. The shown data includes error bars 
taking into account only statistically independent sampleJ^. 

results for the N = 128 system (lower panel), particu¬ 
larly at high depletant reservoir densities (ftp, indicating 
the presence of finite size effects for this system size. 

B. Coexistence curve of the penetrable hard sphere model 

We also tested the capability of our algorithm to equili¬ 
brate hard sphere systems at gas-liquid coexistence, and 
especially near the critical point. We carried out Gibbs 
ensemble simulations of hard spheres in penetrable hard 
sphere depletant^. These types of simulations require 
insertion of the colloid at random positions in the simu¬ 
lation box, which is nearly impossible for high depletant 
fugacities. To overcome this difficulty, we resort to the 
configurational bias scheme discussed in Sec. mE and 
originally introduced in the context of the Gibbs ensem- 
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ble of hard spheres with depleting rods in Ref. [571 For 
every exchange of a colloid between boxes, depletants are 
randomly inserted at the new position, and overlapping 
depletants are attempted to be reinserted in the old box. 
The move is accepted with the probability that accounts 
for the configurational bias weight. 

In Fig. we compare the coexistence curve thus ob¬ 
tained to published data by Vink and HorbactP^I. Those 
authors did not use the Gibbs ensemble, but performed 
direct simulation in the grand-canonical ensemble of the 
colloids and depletants in a single box. Their method is 
advantageous to sample the gas-liquid separation, which 
takes place at intermediate densities (j)c ^ 0.4, because 
it relies exclusively on particle insertion and deletion at 
random positions in the simulation box. Thus, in this 
regime their scheme can be at least as efficient as sin¬ 
gle particle moves, if the particle deletions are combined 
with depletant insertions, and vice versa. However, the 
grand-canonical method is not easily applicable to solid 
phases, for which particle insertion in a crystal lattice is 
nearly impossible. Our method, in contrast, computes 
depletion interactions for single-particle translations and 
rotations. 

As shown in Fig. our data for the total system size 
N = 256, corresponding to the larger of the two system 
sizes studied by Vink and Horbach, generally reproduces 
their data for a depletant-colloid size ratio of g = 0.8, 
at which many-body effects are important. However, we 
see some scatter in our data, which is likely a conse¬ 
quence of surface effects that make it notoriously hard to 
study coexiste nce ne ar the critical point in Gibbs ensem¬ 
ble simulation^SEH Vink and Horbach improved their 
sampling using the umbrella method and thermodynamic 
integration. Overall, however, our data obtained without 
using advanced free energy techniques is in agreement 
with the published data, validating the method. 


V. RESULTS 

A. Aggregation of hemispheres into superlattices 

Equilibrium data of anisotropic particles aggregating 
into crystals with depletants is scarc^^^. Here, we present 
new results on the hierarchical assembly of hemispheres 
into FCC/HCP-cluster phases. Hard hemispheres for 
self-assembly have been the subject of previous investiga¬ 
tion. Marechal and Dijkstra predicted the stability of a 
cluster-FCC (fcc^) phase for hemispheres, but they were 
unable to hnd it in self-assembly simulations of sufficient 
siz^^. Cinacchi presented the phase diagram of hard 
spherical caps, which does not include an fcc^ phas^^. 
Neither study involved depletants. 

We analyze the phase behavior of hemispheres in the 
presence of penetrable hard sphere depletants. Figure 
1^ shows the kinetic phase diagram as a function of de¬ 
pletant reservoir density and colloid density 4>c, for 
a depletant-hemisphere diameter ratio of g = 0.15. Re- 



FIG. 5. Coexistence curve for phase separating hard spheres 
in the presence of penetrable hard sphere depletants. Spheres 
{N = 512) of initial packing fraction (f>c = 0.12 are simu¬ 
lated in the semigrand Gibbs ensemble using at constant nor¬ 
malized depletant reservoir density (j)p = (7r/6)ddep2p using 
implicit depletants (ntriai = 100), and the coexisting colloid 
volume fractions (squares) are obtained fitting the peaks 
of the two-dimensional N — V histogrartP^. Asterisks denote 
data from Ref. [31] measured in the grand-canonical ensemble. 


markably, we observe the formation of the fcc^ and hcp^ 
phases at finite depletant densities > 0.30, and the 
inset shows a snapshot of such a configuration of hemi¬ 
spheres. However, at zero depletant fugacity, which cor¬ 
responds to the case studied previously, we did not ob¬ 
serve any ordered phase, even after 6 x 10® MG sweeps. 
Instead, we find a cluster fluid. In the phase diagram, 
we find close-packed crystals with both HCP and FCC 
stacking, and we suspect the fact that both occur indi¬ 
cates that the free energy difference is smalP^. 

We compare the implicit method against two other 
schemes, an explicit grand-canonical ensemble for the 
depletant^^, and a canonical ensemble with hxed con¬ 
centration of depletants. Figure shows the number of 
hemisphere pairs that have formed after time t. Because 
Monte Carlo simulations do not have a time scale, we 
choose the wall-clock time of the simulation as an ad- 
hoc measure of time. By analyzing bond order, we found 
that the time scale of crystallization corresponds to the 
time when all 512 hemispheres in the simulation box have 
paired up. This event occurs earliest for the implicit 
depletion algorithm. The simulation with explicit grand- 
canonical depletants also orders at a later time. However, 
the simulation with fixed number of depletants does not 
equilibrate into an ordered phase within the wall-clock 
time limit of 48h or 7.8 x 10^ sweeps. Our findings show 
that the implicit algorithm leads to the fastest assembly 
of hemispheres into cluster crystal phases. 
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FIG. 6. Self-assembly of hemispheres into crystalline phases. 
Shown is the kinetic phase diagram for N = 512 hemispheres 
obtained with implicit simulation of depletants as function 
of the depletant reservoir density and the colloid density 
(pc, at depletant-hemisphere diameter ratio q — 0.15. Inset: 
Snapshot of the hcp^ phase found for (pc = 0.575 and = 
0.4. Similar phase diagrams were obtained for q — 0.175 and 
q = 0.125 (not shown). 


B. Diflusivity of discoids with depletants 


Ellipsoids are simple examples of anisotropic particles. 
Recently, discoids have been demonstrated to arrange 
into metastable strand structures at sufficiently high den¬ 
sity of polymeric depletant^. Here, we investigate the 
diffusivity of discoids at depletant densities that do not 
lead to ordering. For Monte Carlo simulations with sin¬ 
gle particle moves, the diffusivity of the colloids in terms 
of mean square displacement per wall clock time is an 
effective measure of the speed of equilibration of the sim¬ 
ulation. In our simulations, we tune the single particle 
step size for translation and rotation so as to yield an 
average acceptance rate of 20%. 

The upper panel of Figure shows the effect of ntriai 
on the diffusivity D of discoids. The colloid particles 
are uniaxial ellipsoids with semi axes a = 6 = 0.5 and 
c = 0.25, the depletants are of radius r = 0.25, and 
the simulations are performed in a dilute system at col¬ 
loid density (pc — 0.01 and depletant reservoir density 
(pp — 0.40, below the coexistence density for metastable 
cluster^^. From the graph, it can be clearly seen that 
using configurational bias moves with a modest value of 
’T-triai ^ 10 speeds up the equilibration by almost three 
orders of magnitude compared to not using configura¬ 
tional bias moves. The effect is dramatic and similar in 
magnitude between running the simulation on the CPU 
vs. the GPU. At peak diffusivity, there is a slight ad¬ 
vantage to using the GPU, compared to CPU socket per¬ 
formance. For higher values of ntriai, tbe performance 
drops off slowly, as a result of the increased computa¬ 
tional effort to carry out the depletant reinsertions, while 
the effect of increasing the step size due to a higher ac¬ 


FIG. 7. Aggregation kinetics of hemispheres. Shown is 
the nnmber of spheres formed after simulation time t (in 
hours), for a simulation with implicit depletants (diamonds), 
explicit grand-canonical depletants (squares) and canonical 
depletants (circles). Simulations where performed at colloid 
volume fraction pc = 0.575 and depletant reservoir density 
Pp = 0.40 for a depletant-colloid diameter ratio of g = 0.175, 
on eight cores of an Intel Xeon E5-2680 processor with spatial 
domain decomposition via MPI (single precision). A sphere is 
defined as two hemispheres with their face centers being closer 
than 0.2d apart, where d is the diameter of the (hemi-)sphere. 
In the canonical case, the constant number Np = 4884 of 
explicit depletant particles has been chosen to be the aver¬ 
age number of depletants in the free volume of the grand- 
canonical simulations, after phase transformation. 


ceptance ratio is weaker. We note that we carried out 
simulations with finite values of ntriai at higher colloid 
densities as well (data not shown) and found the effect 
to be less pronounced at these densities. 


We further measure the performance at different col¬ 
loid densities pc between the dilute regime and the regime 
of a dense liquid, for the same parameters as above, with 
?T-triai = 0 (Fig. i lower panel). For simulations with 
implicit depletants, either using the CPU or the GPU, 
the performance depends only slightly on the colloid vol¬ 
ume fraction, directly confirming the beneficial effect of 
implicit calculation of the interaction in the dilute sys¬ 
tem, where the number of depletants would be very high 
with an explicit treatment. Indeed, the performance of 
the explicit depletant simulations in the grand-canonical 
ensemble drops noticeably when going from pc = 0.50 
towards lower densities, and the system becomes practi¬ 
cally impossible to equilibrate when pc < 0.30. Looking 
at GPU vs. CPU performance, we note that GPUs are 
advantageous for very dilute systems, but do not pro¬ 
vide better performance when the system is dense in col¬ 
loids. This is because the checkerboard parallelization 
scheme implemented for performing the colloid moves on 
the GPU (Sec. HIE and Ref. |33) requires a large simu¬ 
lation box to operate efficiently. 
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FIG. 8. Diffusivity of discoids in penetrable hard sphere de- 
pletants. Upper panel: Diffnsion coefficient vs. the num¬ 
ber ntriai of configurational bias swaps, for a simulation of 
N = 500 discoids on 12 CPU cores (Intel Xeon E5-2680v3) ns- 
ing MPI (squares) and a single NVIDIA K20X GPU (circles). 
For simnlation parameters, see main text. The diffnsivity 
is obtained from fitting the linear mean square displacement 
MSD as function of the wall-clock time t (in seconds). Lower 
panel: Diffnsion coefficient vs. colloid density (pc (utriai = 0), 
for a simulation on a NVIDIA K80 GPU (sqnares), on 8 cores 
of an Intel Xeon E5-2680v3 (circles), and for a simulation 
of explicit depletants in the grand-canonical ensemble (dia¬ 
monds), on the same hardware. 


VI. CONCLUSION 


We see applications for our method in the simulation 
of anisotropic colloid phase behavior. Even without de¬ 
pletants, polyhedra have been shown to order into a 
multitude of different structured. With dep letion in- 
teractions, additional phases can be stabilize d^ l ^^ l ^^ l ^'^ . 
The algorithm can also be used to study the aggrega¬ 
tion of entropically patchy colloids into colloidal polymer 
chains, held together by strong depletion bonda^^. In 
this context, it would be interesting to study solutions 
as well as melts of such colloidal polymers. An inter¬ 
esting open question concerns whether depletant entropy 
can stabilize not only close-packed but also open ordered 
structured^. In protein crystallization, depletant poly¬ 
mers are commonly used as precipitants. An important 
limitation of our algorithm is that it treats only non¬ 
interacting depletants, and the validity of that approx¬ 
imation remains to be investigated for specific systems. 
In contrast to enthalpically patchy models, our algorithm 
does not require implementation of shape-specific attrac¬ 
tive patches to study aggregation of colloids, and the al¬ 
gorithm is therefore highly robust and generic. 
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We have presented an efficient algorithm to implicitly 
simulate depletion interactions between anisotropic col¬ 
loids. The algorithm is implemented on parallel multi¬ 
core processors and graphics processing units. Com¬ 
bined with a parallel Monte Carlo schem^^mH^ 
gorithm offers a way to tackle large scale simulations of 
hard shapes with depletants. The scheme may be readily 
generalized to soft interactions between the colloid and 
the depletant, such as the Hertz potentiaP^. We stress 
that even though the algorithm is parallel, already its se¬ 
rial implementation offers significant speed-ups over algo¬ 
rithms that do not use cluster moves, for dilute systems of 
colloids, because only depletants in the neighborhood of 
every particle are considered. Nevertheless, the method 
works perfectly well for the fluid-solid transition. 


Appendix A: GPU implementation 

In the GPU implementatiom we perform the colloid 
trial moves in the active cells^ and the depletant inser¬ 
tions in different kernels. To insert depletants, we draw 
a random number of depletants for every moved colloid, 
as described in Sec. |HIB| We use a one-to-one mapping 
between depletants and thread groups of size n < nmax- 
Here, rimax = 32 is the maximum number of threads that 
can perform overlap checks synchronously, and we tune 
n at run-time. When any thread detects an overlap be¬ 
tween the depletant and any particles in the old configu¬ 
ration, the depletant is ignored. In the other case, if the 
depletant overlaps with the moved colloid, that colloid 
move is flagged for rejection. 
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When the configurational bias scheme is used (ntuai > 
0), a second kernel with a similar thread mapping is 
launched, however, depletants are assigned to whole 
thread blocks of size s < 1024, which is an auto-tuned 
parameter, so that the bias weights of difi’erent reinser¬ 
tions belonging to the same depletant can be summed in 
shared memory. 
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